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We present a concise, but systematic, review of the ergodicity issue in strongly 
correlated systems. After giving a brief historical overview, we analyze the issue 
within the Green's function formalism by means of the equations of motion ap- 
proach. By means of this analysis, we are able to individuate the primary source 
of non-ergodic dynamics for a generic operator and also to give a recipe to 
compute unknown quantities characterizing such a behavior within the Com- 
posite Operator Method. Finally, we present examples of non-trivial strongly 
correlated systems where it is possible to find a non-ergodic behavior. 
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1. Historical Overview 

The issue of ergodicity in condensed matter physics is well known since fifties [ 
[1]. Given two operators A and B, describing physical quantities (e.g., charge, spin, 
pair densities or currents), one can study the physical response of a system described 
by a certain Hamiltonian H through the generalized susceptibility 

XAB = lim^ (1.1) 

h-^o on 



where h is an external field entering the Hamiltonian of the system under study 
in a coupling term of the type —hB and (■ ■ ■ ) stands for the statistical average in 
some ensemble over the perturbed system. Kubo [ [1] immediately noticed that the 
static isolated susceptibility x^(0), defined for an isolated system perturbed by an 
external field turned on adiabatically, and the isothermal susceptibility defined 
for a system in thermal equilibrium in presence of a time- independent external field, 
can be generally different. In particular, Falk [[2] has shown that the static isolated 
susceptibility is just a lower bound for the isothermal one x^(0) < ■ 
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We recall that the static isolated susceptibility x^(0)) within the linear response 
theory, is defined through the related retarded Green's function [[1] 

X^M = -T \vB{h - mu), Bit,)])^] (1.2) 

where (■ ■ ■)o stands for the statistical average in the microcanonical ensemble (i.e., 
fixing energy) on the unperturbed system and JF for the Fourier transform. 
On the other hand, the isothermal susceptibility can be computed as 

[\Ai-i\)B)od\-P{A)o{B)o (1.3) 

In fact, starting from the expression of the thermal average in the canonical ensemble 
(i.e., fixing temperature) (A) 

(A) = iTr(Ae-'3(^"^^)) (1.4) 

where /3 = |; and Z = Ti (^^-PH-hb'^ ^ ^^^^ expand ^-PiH-hB) powers of h and 
get 

^-m-hs) ^ ^-pH (^ + hj^^ dAe^^Ee-^^ + 0(/i2)^ (1.5) 

Substituting this expansion into (11.41) and retaining only the first order term in 
we get at the numerator 

Tr(Ae-^(^-'^^)) ^Zo(A)o + /iZo / d\{AB{i\))o (1.6) 

Jo 

and at the denominator 

Z^Z^{l + hp{B)^) (1.7) 

and by the ratio 

{A) ^ (A)o + h f d\{AB{i\))^ - hp{A)o{B)o (1.8) 



where Zq denotes the partition function of the unperturbed system. Taking the 
derivative after (11.11) and exploiting the cyclic property of the trace, we obtain the 
isothermal susceptibility as in (11.30 . 

Now, if we rewrite both expressions by means of the general formulas for the 
retarded Green's functions and the correlation functions given in the companion 
article [[3] (see Section 3), present in this same issue, we get 



and 



= E rAB(k) ^\-Y^ - muB), (1.10) 
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where [ E], H] V is proportional to the volume of the system, the sum over / ranges 
over the number of fields in the chosen basis, a'-'-' are the spectral density functions, 
a;^^) are the poles of the propagator. Tab is an unknown function appearing in case 
of poles with zero value. 

We can immediately see that the two susceptibilities differ for the following 
expression 

- X\0) = pIpY1 r^^(k) - P{A)o{B)o (1.11) 

k 

Now, one can check that rewriting the expression limt^^{AB{t)) by means of 
the general formula for the correlation functions given in the companion article [ E] 
(see Section 3), present in this same issue, we just get 

}im{AB{t)) = ^Y.^ABik) (1.12) 

k 

and, accordingly, the difference at the r.h.s. of (11.111) is just what enters Khintchin's 
theorem [[5]: a dynamics is ergodic (i.e., phase space equilibrium averages are equal 
to ensemble microcanical averages, which are much easier to compute) 

/"OO 

{AB) = / dtAB{t) (1.13) 
Jo 

if an only if 

\mi{AB{t)) = {A){B) (1.14) 

t— >oo 

In other words a dynamic is ergodic if correlations attenuate in time. In particular, 
for B = A, the dynamics of A is ergodic if, during its time evolution, it has non-zero 
matrix elements only between states within a zero-volume region of the phase space 
of the system [ E] . 

It is clear now the link between ergodicity and response theory: the two defini- 
tions of susceptibility differ when the dynamics of the system is not ergodic. Two 
little, but important, notices: finite systems are not ergodic by definition, just be- 
cause of the inequivalence of the ensembles; non-ergodicity at zero temperature is 
just the result of a degeneracy in the ground state. 

Several years later it was shown [[71[8] that the difference between the two defini- 
tions of susceptibility is related to the zero-frequency anomaly exhibited by bosonic 
correlation functions: the presence of undetermined constants in the bosonic corre- 
lation functions. This is exactly what the relations derived above and the results of 
the companion article [|3] predict establishing a definite link between the ergodicity 
of the dynamics and the Green's function formalism. It was first put in evidence in [ 
E] and then studied by many other authors [li0l[IIl[71[8l|6l[l2l[T3l[Tl[l5l[l6]. There 
is a general believe that this problem is of academic interest and in the last years no 
much attention has been dedicated to it. The main reason is that the response func- 
tions, the experimentally observed quantities, are given by retarded bosonic Green's 
function which formally do not depend on the such undetermined constants, which 
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are, therefore, considered of no physical interest. The general attitude [ [H, is 
to believe that in macroscopic real systems at equilibrium at a temperature T, the 
fluctuations are very small and the interaction between the system and the reservoir 
would introduce an irreversible relaxation and decouple the correlation functions. 
Then, as suggested in [ (TU] , these constants should be always determined by requir- 
ing the ergodicity. This procedure is some how an artifice and may lead to serious 
problems because it might break the internal self-consistency of the entire formu- 
lation. As remarked in [[10], the zero-frequency anomaly is a manifestation of the 
difficulty in extracting irreversible behavior from the statistical mechanics. This is 
true, but as long as we use the scheme of statistical mechanics we must be careful 
in doing self-consistent calculations. Breaking the self-consistency might bring to 
serious errors. 

According to the well-known relations existing between casual (C), retarded {R) 
Green's functions and correlation functions 

3f?[G^(k, u)] = 3?[G^(k, to)] (1.15) 



^^[^"(k,^;)] = tanh ( ^ ) 55[G^(k,cu)] (1.16) 



C(k,cu) 



1 + tanh I 



»[G^(k,..)] (1.17) 



the zero-frequency excitations do not contribute explicitly to the imaginary part 
of the retarded Green's functions and, consequently, T does not explicitly appear 
in the expressions of susceptibilities. At any rate, susceptibilities retain an implicit 
dependence on F through the matrix elements. Then, the right procedure to compute 
both correlation functions and susceptibilities is clearly the one that starts from the 
causal Green's function, which is the only Green's function that explicitly depends on 
r. It is worth noticing that the value of F dramatically affects the values of directly 
measurable quantities (e.g., compressibihty, specific heat, magnetic susceptibility, 
. . . ) through the values of correlation functions and susceptibilities. According to 
this, whenever it is possible, F should be exactly calculated case by case. 

If we do not have access to the complete set of eigenstates and eigenvalues of 
the system, which is the rule in the most interesting cases, we have to compute cor- 
relation functions and susceptibilities within some, often approximated, analytical 
framework. Now, since no analytical tool can easily determine F (e.g., the equations 
of motion cannot be used to fix F as it is constant in time), one usually assumes the 
ergodicity of the dynamics of ip and simply substitutes F by its ergodic value (i.e., 
by the r.h.s. of ffLTi]) ): 

r-^{i,i) = mi)){^^m- (1-18) 

Unfortunately, this procedure cannot be justified a priori (i.e., without computing 
F through its definition (12.581) ) by absolutely no means. The existence of just one 
integral of motion and, more generally, of any operator that has a diagonal part 
with respect to the Hamiltonian [ [6] (i.e. by any operator that has a diagonal en- 
tries whenever written in the basis of eigenstates of the Hamiltonian) divides the 



4 



Ergodicity in Strongly Correlated Systems 



phase space into separate subspaces not connected by the dynamics. This latter, 
in turn, becomes non ergodic: time averages give different results with respect to 
ensemble averages. This latter consideration also clarifies why the ergodic nature of 
the dynamics of an operator mainly depends on the Hamiltonian it is subject to. 

It is really remarkable that F is directly related to relevant measurable quantities 
such as compressibility and specific heat trough the dissipation-fluctuation theorem. 
For instance, we recall the formula that relates the compressibility to the total 
particle number fluctuations 

(1.19) 

where N is the total particle number operator, is its average and V is proportional 
to the volume of the system. We see that a compressibility different from zero requires 
the non-ergodicity of the system with respect to total particle number operator. 
According to this, in the case of infinite systems too the correct determination of F 
cannot be considered as an irrelevant issue (e.g., fll.lQp holds in the thermodynamic 
limit too). 

In the next section, we provide some examples of violation of the ergodic con- 
dition fll.l4p . It is necessary pointing out, in order to avoid any possible confusion 
to the reader, that we are using full operators and not fluctuation ones (i.e., we use 
operators not diminished of their average value, in contrast with what it is usually 
done for the bosonic excitations like spin, charge and pair). According to this, the 
F can be different from zero (i.e., be equal to the squared average of the operator), 
and still indicate an ergodic dynamics for the operator. 



V 

Ar2 I 



2. Examples 

2.1. Two-site Hubbard model 

The two-site Hubbard model is described by the following Hamiltonian 

H = Y1 (^'^- - ^^^^') c^(^) + ^E^t(0 n^it) (2.1) 

ij i 

where the summation range only over two sites at distance a from each other and 
the rest of notation is standard [H]. The hopping matrix tij is defined by 

= -2t a,, a,, = ^Y1 ^''^"'^ «(^) (2.2) 

k 

where a{k) = cos{ka) and k = 0, n/a. 

We now proceed to study the system by means of the equation of motion ap- 
proach and the Green's function formalism [ |T7] . A complete set of fermionic eigen- 
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operators of if is the following one 



r]{i) 



(2.3) 



where 



m 

r]{i) 



[1 - n{^)] c{i) 



nn) cn) 



(2.4a) 
(2.4b) 

(2.4c) 
(2.4d) 



We define ip'^ 



J2j (^ij ^(j) use the spinorial notation for the field operators. 



n/i(z) = (i) a^c{i) is the charge (/x = 0) and spin (/i = 1, 2, 3) operator; greek 
(e.g., /i, u) and latin (e.g., a, b, k) indices take integer values from to 3 and from 
1 to 3, respectively; sum over repeated indices, if not explicitly otherwise stated, is 



understood; cr„ 



[1, a) and cr^ 



-l,a): a are the Pauli matrices. In momentum 



space the field satisfies the equation of motion 
where the energy matrix e(fc) has the expression 



(2.5) 



e{k) 



/ -/i - 2t(x{k) -2ta{k) 

f/-/i 

At 

\ 2t 



-2t 
2t 



-2t 
2t 



\ 



-fi + 2ta{k) Ata{k) 
2ta{k) U - fi ) 



(2.6) 



Straightforward calculations, according to the scheme traced in [ [T7] , show that 
two correlators 



(2.7) 
f2.8) 



appear in the normalization matrix J(k) = T ({'?A(i, t), '?/'^(j, t^\)- Then, the Green's 
functions depend on three parameters: /i, A and p. The correlator A can be expressed 
in terms of the fermionic correlation function C{i,j) = 'ip'^ {j)y, the chemical 
potential fi can be related to the particle density by means of the relation n = 
2 [1 — Cu{i,i) — C22{'i,i)]- The parameter p cannot be calculated in the fermionic 
sector; it is expressed in terms of correlation functions of the bosonic fields nf^{i) and 
ci{i)ci{i). According to this, the determination of the fermionic Green's functions 
requires the parallel study of bosonic Green's functions. 
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After quite cumbersome calculations, it is possible to see [[T7] that a complete 
set of bosonic eigenoperators of H in the spin-charge channel is given by 



where 



with the definitions: 



Bf{^=d,{^-dl{{)+dl{i)-dl-{^ 
Bi'\z) = d,{i)-d;{z)~dl{z) + dl-{i) 

B^''\^) = U^)-f^{^ + fl{^-fl''i^) 



The field B'^^\i) satisfies the equation of motion 

i^B^^\k,t) = K{k) B^^^k^t) 
where the energy matrix K,{k) has the expression 



K{k) 








-2t 











\ 




-4t[l -a(fc)] 





u 






















u 


2t 













u 








2t 










8t 











v 











8t 





/ 



The energy spectra are given by 



ui{k) 


= -2t^2 [1 


-a{k)] 


LU2{k) 


= 2t^2 [1 - 


a{k)] 


UJ3{k) 


= -f/ - 




U4{k) 


= -4:JU 




u^{k) 


= ^Ju 






= U + AJu 





(2.9) 



(2.10) 
(2.11) 
(2.12) 
(2.13) 
(2.14) 
(2.15) 

(2.16) 
(2.17) 

(2.18) 



(2.19) 



(2.20) 



(2.21) 

(2.22) 
(2.23) 
(2.24) 
(2.25) 
(2.26) 
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where 



Ji 



u 



(2.27) 



Straightforward calculations show that the correlation function has the expres- 



sion 



C(^)(z,j) = (5(^n0 5^^^^(j)> 
1 ^ 

k n=l 



k(i-j)-i(jJn{k){ti-tj) 



1 + tanh 



/("''^)(fc) (2.2^ 



where 



j(n,M)(0) 
j{n,ii) 



TT 



for n = 3, 4, 5, 6 

coth^^a(-'^)(vr) 
2 ^ ' 



Owing to the fact that zero-energy modes appear for n = 1, 
Eq. fl2.2ip ]. r appear in the correlation functions 



1 ^ 

r(^Ho) = -5^/("''^)(o) 



(2.29a) 
(2.29b) 

2 and k = Q [cfr. 

(2.30) 



n=l 



One might think, as is often done in the literature, to fix this constant by its ergodic 
value. However, this is not correct as we are in a finite system in the grandcanonical 
ensemble and the ergodicity condition does not hold. For the moment, we can state 
that this constant remains undetermined. 

The spectral density functions depends on a set of parameters which come from 
the calculation of the normalization matrix I^^\k) = J-" (^[B^^\i,t), 5*-'^-'^(j, t)] ). In 
particular, for the (l,l)-component the following parameters appear: 



C- 

d-- 



12 



U c 



c\{i)c\{i) 



(2.31a) 
(2.31b) 

(2.31c) 

(2.31d) 



The parameters and C^2 related to the fermionic correlation function C{i,j) = 
{ipi^i) 'ip^ij))- The parameter Xs can be expressed in terms of the bosonic correlation 
function C^^\i,j) = {^B^^\i) B^^'>\j)') . In order to use the standard procedure of 
self-consistency, we need to calculate the parameter d. For this purpose we should 
open both the pair channel and a double occupancy-charge channel (i.e., we will 
need the static correlation function {n^{i)ni{i)n°'{i))). The corresponding calcula- 
tions are reported in Ref. [[T7| where is shown that these two channels do not carry 
any new unknown F. The self-consistence scheme closes; by considering the four 
channels (i.e., fermionic, spin-charge, pair and double occupancy-charge) we can set 
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up a system of coupled self-consistent equations for all the parameters. However, 
r(/^)(^0) has not been determined yet: we have not definitely fixed the representation 
of the Green's functions. 

In conclusion, the standard procedure of self-consistency is very involved and is 
not able to give a final answer because of the problem of fixing the F. We will now 
approach the problem by taking a different point of view. The proper representation 
of the Green's functions must satisfy the condition that all the microscopic laws, ex- 
pressed as relations among operators must hold also at macroscopic level as relations 
among matrix elements. For instance, let us consider the fermionic channel. We have 
seen that there exists the parameter p, not explicitly related to the fermionic prop- 
agator, that can be determined by opening other channels. However, we know that 
at the end of the calculations, if the representation is the right one, the parameter p 
must take a value such that the symmetries are conserved. By imposing the algebra 
constraints (??) and by recalling the expression for A we get three equations 

n = 2(1 - Cn - C22) (2.32a) 
A = Cfi - C^, (2.32b) 
Cu = (2.32c) 

This set of coupled self-consistent equations will allow us to completely determine 
the fermionic Green's functions. Calculations show [ |T7j that this way of fixing the 
representation is the right one: all the symmetry relations are satisfied and all the 
results exactly agree with those obtained by means of Exact Diagonalization. We do 
not have to open the bosonic channels; the fermionic one is self-contained. 

Next, let us consider the spin-charge Green's functions. In the spin-charge sector 
we have the parameters C", 0^2, ? d and two F 

^o = iE/n°^(0) (2.33) 

i=l 

h = \i2fu^(^) fc = 1,2,3 (2.34) 

i=l 

Since we are in absence of an external applied magnetic field, bk takes the same 
values for any value of k. 

The parameter and Cfg are known, since the fermionic correlation functions 
have been computed. The parameters Xs d, can be computed by means of the 
equations 

d=^(n^(z)n^(0>-p (2.35) 
= (n(0-n"(^)) (2.36) 

The F are fixed by the algebra constraints 

Cit\z,z) = {n,{z)n,{z)) (2.37) 
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Figure 1. (left) 6o and are plotted as functions of n for [7 = 4 and T = and 
1. U and T are expressed in units of t. (right) 6o and bk are plotted as functions 
of U for T = 0.01 and n = 0.6, 0.8, and 0.9. U and T are expressed in units of t. 



By recalling fl^:^ and (ET^ we have 



1 ^ 

E 



i=l 



1 + coth 



(3Ui{'K) 



a 



11 ' 



(2.38) 



with 



n + 2D for jj, 
n — 2D for /i 



1,2,3 



(2.39) 



D = {n^{i) ni{i)) is the double occupancy and can be calculated by means of the 
fermionic correlation functions D = n — 1 + Cu. Eqs. f l2.35p and fl2.38p constitute a 
set of coupled self-consistent equations which will determine completely the Green's 
function in the spin-charge channel. Calculations show that this way of fixing the 
representation is the right one: all the symmetry relations are satisfied and all the 
results exactly agree with those obtained by means of Exact Diagonalization. 

bo and bk are plotted as functions of n and U in Fig. [T] for various temperatures. 
It is worth noting that they assume their ergodic values (i.e. and 0, respectively) 
only in some regions of the parameter space: (at zero temperature) at n = 1 (both 
bo and bk) and at n = 0.5 (5o only). In these regions, the grand-canonical ensemble 
is equivalent to the microcanonical one and the underlying ergodicity of the charge 
and spin dynamics emerges. 

It is worth noting that bo is directly related to the compressibility by means of 
the following relation [ [17] 



According to this, if we erroneously set the value of bo to the ergodic one (i.e., n^) 
we would get a constant zero compressibility. 
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2.2. Tight-binding model 



A narrow-band Bloch system in presence of an external magnetic field is de- 
scribed by the following Hamiltonian 



H = J2{Ui-ii ct (z) c{j) -hY, nsiz) 



(2.41) 



where 713(2) is the third component of the spin density operator and h is the intensity 
of the external magnetic field. The indices i and j run on an infinite rf-dimcnsional lat- 
tice. Straightforward calculations show that the causal Green's function G^''(i, j) = 
(T [n^(i) nfj,{j)]) and the correlation function C'^^^ (i, j) = (n^(i) n^{j)) of the charge- 
spin operator n^{i) — c^{i) a^c{i) have the following expressions 



C^'^Hk, c^) = (27r)'^+^ a"^ 5^'^\k) 5{u) V^^'^ + 



1 

1 -I- tanh -— 



(2.42) 

a[Q(''H k,a;)] (2.43) 



where ^'^^{k) is the d-dimensional Dirac delta function. (5^''^(k, a;) comes from the 
proper fermionic loop and is the Fourier transform of 



g('^)(z,j) = Tr[a^Gc(i,J>MGc(jV 



(2.44) 



Here Gcii-,]) = (T [c{i) (j)]) is the causal fermionic function and has the 
expression 



Gc(k,a;) = J]- 



a 



in) 



n=l 



g-/3E„(k) 



1 



g-/3i?„(k) 



u - En(k) +16 uj- En{k) - iS 



with 



(2.45) 



where 



^i(k) 



<7« = 



—fj, — 2dta(k) — h 
-ljL-2dta(\s.) + h 
1 0\ ^(2)^fO 

vol 



1 

aCk) = - cos{ki 
1=1 



is fixed by the algebra constraints (??) which requires 



^('^) = (n^(^)r^^(^))- 



(27r) 



d kduj 



1 -|- tanh ■ 



(2.46) 
(2.47) 

(2.48) 



(2.49) 



»[Q(^)(k,cu)] (2.50) 
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The loop (5^'^'*(k, tu) can be calculated by means of fl2.45p . Calculations show 



d+l 



d'^k duj 



1 13^' 
1 + tanh 



= (n)-K)'-K)' /or ^ = 0,3 (2.51) 
= {n)-2{n^{i)){n^{^) for fi = l,2 (2.52) 

By recalling the algebra constraints (12.391) . Eq. (12.501) gives for the F 

r(°) = {nf (2.53) 
r^^'^) = (2.54) 
r(^) = {n^f (2.55) 

in accordance with the ergodic nature of the spin and charge dynamics in this system. 

It is worth noting that the compressibility of this system can be computed by 
means of the general formula ( I1.19P that holds in the thermodynamic limit too and 
gives 



1 p a" 



E 



d'^k 



C„(k) 



(2.56) 



{nf 2 2(27r)'^ 

where C„(k) = cosh^ (^^Mlj. We can see that an ergodic charge dynamics can 

lead to a non-ergodic value of the F relatively to the total number operator, which 
is an integral of motion. Also in the infinite systems the decoupling inspired by the 
requirement of ergodicity cannot always be applied. 



2.3. Heisenberg chain 

We will now study [HH [19] the ergodicity of the dynamics of the operator S?, 
the z-component of the spin at site i, in the ID anisotropic extended Heisenberg 
model described by the following Hamiltonian: 

H = -J^J2 S-SUi + J± + ^i'^i+i) + ^' E (2.57) 

i i i 

where S"?, Sf and S? are the x, y and z components of the spin-1/2 at site i, 
respectively. The model (12.571) is taken on a linear chain with periodic boundary 
conditions. We take the interaction term parameterized with ferromagnetic ( > 
0) and the next-nearest-neighbor interaction term, which is parameterized with J', 
isotropic. In order to frustrate ferromagnetism, we have considered only the case with 
J' > 0, that is, with an antiferromagnetic coupling between next-nearest neighbors. 
According to this, only chains with even number of sites have been studied in order 
to avoid topological frustration that would be absent in the thermodynamic limit. 
Since it is possible to exactly map all results obtained for > to those for < 
by means of a simple canonical transformation, we have limited our study only to 
positive values of J±. 
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We have numerically diagonalized the Hamiltonian (12.571) for chains of size L 
ranging between 6 and 18 by means of Exact Diagonalization (ED) (divide and con- 
quer algorithm) and for chains of size L ranging between 20 and 26 by means of 
Lanczos Diagonalization (LD). We have systematically taken into account transla- 
tional symmetry and classified the eigenstates by the average value of S*^ = 5?, 
which is a conserved quantity. Whenever we have used ED, all eigenvalues and eigen- 
vectors of (12.571) have been calculated up to machine precision and, therefore, we have 
been able to determine the exact dynamics of the system for all temperatures. On 
the contrary, when we have used LD, we have been limited to the zero-temperature 
case since only the ground state can be considered exact in LD. 

In this case, we have the opportunity to exactly compute F in terms of the exact 
eigenvalues and eigenstates \n) of the system. As a matter of fact, it read as [0] 

r = ^ ^ e''^^-{n\Sl\ni){ni\Sl\n) (2.58) 

EfL—Em 

As already discussed above, the dynamics of an operator (e.g., 5?) is ergodic when- 
ever (11.141) is satisfied, or equivalently, (12.581) is equal to its ergodic value: 

^erg ^ ^^z^2 = i_ ^ e-^(^"+^'") (n 1 5f | n) (m | Sflm) . (2.59) 

n,m 

The dynamics of a finite system is hardly ergodic, since (I2.58P and (12.591) unlikely 
coincide. In the thermodynamic limit, the sums in (12.581) and (12.591) become series 
and no conclusion can be drawn a priori. Since we have diagonalized the Hamiltonian 
(12.571) numerically (i.e., only for finite systems) and since L — cx) is the most 
interesting case, we have analyzed our results through finite-size scaling in order to 
speculate on the properties of the bulk system. 

If the ground state of fl237l) is A^-fold degenerate then, at T = 0, fl238l) and fl239l) 
read as follows: 



n,m=l 




respectively. 

Thanks to the translational invariance enjoined by the system (S?) is independent 
of i and proportional to the ^r-component of the total spin operator average (S^^^). It 
is easy to show that, even if there is a finite magnetic moment per site in any of such 
N degenerate ground states, (S'f ) at T = is always zero in absence of magnetic 
field. Indeed, if a ground state with non-zero (S'f) = M exists, also another ground 
state with (S'f) = —M exists. Thus, at zero temperature, V^^^ is always zero and 
the only quantity of interest is F. A finite value of this latter implies non-ergodicity. 
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Figure 2. Zero-temperature ergodicity phase diagram in the J' — J± plane. Due 
to the symmetry of the Hamiltonian only the upper half is shown (see in the 
text). Only two ergodic phases (E-I and E-II) have been found in the reported 
parameter space. The others are either non-ergodic (NE-I and NE-II) or impos- 
sible to conclusively analyze (W). The latter phase might shrink to a transition 
line in the bulk limit. 



Obviously, if = 1 then both values coincide. Therefore, a non-ergodic phase 
corresponds to degenerate ground states with finite magnetization. 

In the studied range of coupling constants (see Fig. [21) we have found two non- 
ergodic phases (NE-I and NE-II), two ergodic ones (E-I and E-II) and a weird phase 
(W). Our computational facilities limit the range of chain sizes that we can analyze 
such that we could not establish, by means of finite-size scaling, whether the weird 
phase (W) is ergodic or not. In the non-ergodic phases (NE-I and NI5-II), we were 
able not only to perform the finite-size scaling, but also to write down an analytic 
expression for F as a function of the chain size L. The weird phase (W) has exhibited 
a strong dependence of the ground state upon the particular values of the couplings. 
On the contrary, the other phases exhibit ground states that are independent of the 
particular values of the coupling constants. 

In the standard Heisenberg model (J' = and Jj_ = Jz) at T = the dynamics 
is non-ergodic for ferromagnetic coupling ( J± = < 0) as the system has a L -|- 1 
degenerate ground state 

F = — + — . (2.60) 

12 6L ^ ' 

It is clear from fl2.60p that F remains non-ergodic also in the thermodynamic limit. 
This point (J' = and J± = Jz) becomes a hne in our phase diagram and is denoted 
as NE-I (see Fig. [2]). In fact, the next-nearest-neighbor interaction J' may frustrate 
(J' > 0) or favor (J' < 0) the ferromagnetism. In the latter case, the ground state 
remains unchanged for any value of J' < 0. Therefore, we expect the line denoting 
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Figure 3. Finite-size scaling in the case of T = for different points in the phase 
diagram of Fig. [2l Symbols on panel a): + corresponds to (NE-II), ▲ corresponds 
to (NE-I), ■ corresponds to (E-I) and Q to (E-II) regions of Fig. [21 respectively. 
On panel b) different examples from (W) region are shown. Hamiltonian couplings 
are shown in the legend. All energies are expressed in units of J^. 



the phase NE-I to extend also to negative J' . If, on the contrary, J' is positive and 
large enough to frustrate the system in such a way that the ground state loses its 
ferromagnetic character, the ergodicity is restored. This occurs at a finite critical 
J' ~ O.25J2. For values of J' larger than the critical one, we find a non-degenerate 
ground state with {S^^^) = 0. 

If J_|_ 7^ Jz the rotational invariance is broken so that states with the same (S|,J, 
but different (S^^^), are not degenerate anymore. In the non-ergodic region (NE-II) 
of the phase diagram (see Fig. [2]), the ground state is just doubly degenerate (not 
L + 1 degenerate as in (NE-I)): one ground state corresponds to a configuration with 
all spins up and the other to a configuration with all spins down. Hence, the value 
of F in this phase is 1/4 and does not depend neither on the Hamiltonian couplings 
nor on the number of sites in the chain. It is clear that also this phase extends to 
negative values of J'. This kind of ground state stands the frustration introduced 
by next-nearest-neighbor interaction up to J' ~ 0.3 Jz (see Fig. [2]). 

The ergodic region (E-I) of the phase diagram (see Fig. [2]) has F = for all sizes 
of the system and values of the couplings: the unique ground state belongs to the 
sector with (S^^f.) = 0. On the contrary, the other ergodic phase (E-II) has non-zero 
values of F for values of L not multiples of four. The ground state in this phase has 
average total spin equal to one and, therefore, F = 1/L^. We obviously conclude 
that (E-II) phase is ergodic in the thermodynamic limit. 

The values of F in these four phases (NE-I, NE-II, E-I and E-II) exhibit per- 
fect finite-size scaling as shown in Fig. [3^). This has allowed us to make definite 
statements also in the thermodynamic limit. 

The weird phase (W) (see Fig. [2]) is characterized by a quite strong size depen- 
dence, as shown in Fig. ^) where a tentative finite-size scaling of F in the different 
points of the phase is presented. This region manifests a diverging finite-size scaling 
within the range of sizes we were able to handle. In this case, the behavior of F 
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as a function of L strongly depends on the particular choice of the Hamiltonian 
couplings and is highly non monotonous when increasing L, according to the strong 
dependence on L of (S^^^) in the ground state. In this critical region the eigenval- 
ues of (12.571) present many level crossings, which means that the maximum value 
of L we were able to reach {L^ax = 26) is not large enough to perform a sensible 
finite-size scaling analysis. However, we expect that this phase becomes ergodic in 
the thermodynamic limit, although still different from the ergodic phases E-I and 
E-II. 

We can summarize our findings in the thermodynamic limit at zero temperature 
as follows: 

if J± = ±J, and J' < 0.25J, 



j_ 

12 
1 
4 

??? 



if\J±\< J, and J' < 0.3 J, 

in the weird phase (W) (see Fig. [2]) 

otherwise 



(2.61) 



3. Conclusions 



In conclusion, we have analyzed the issue of ergodicity, after a brief historical 
overview, within the Green's function formalism by means of the equations of mo- 
tion approach. We have individuated the primary source of non-ergodic dynamics 
for a generic operator in the appearance of zero-frequency anomaly in its correla- 
tion functions and given a recipe to compute the unknown quantities characterizing 
such a behavior within the Composite Operator Method. Finally, we have presented 
examples of non-trivial strongly correlated systems where it is possible to examine 
a non-ergodic behavior: two-site Hubbard model, tight-binding model, Heisenberg 
chain. 
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